Theoretical Study of the Structure and Binding Energies of Dimers of Zn(II)-Porphyrin Derivatives

Zinc-complexed porphyrin and chlorophyll derivatives form functional aggregates with remarkable photophysical and optoelectronic properties. Understanding the type and strength of intermolecular interactions between these molecules is essential for designing new materials with desired morphology and functionality. The dimer interactions of a molecular set composed of porphyrin derivatives obtained by substitutional changes starting from free-base porphyrin is studied. It is found that the B97M-rV/def2-TZVP level of theory provides a good compromise between the accuracy and cost to get the dimer geometries and interaction energies (IEs). The neglect of the relaxation energy due to the change in the monomer configurations upon complex formation causes a more significant error than the basis set superposition error. The metal complexation increases the binding energy by about −6 to −8 kcal/mol, and the introduction of keto and hydroxy groups further stabilizes the dimers by about −20 kcal/mol. Although the saturation of one of the pyrrol double bonds does not change the IE, the addition of R groups increases it.


■ INTRODUCTION
Metalloporphyrins are ubiquitous in light-harvesting organisms. 1 Because of their extraordinary photophysical and optoelectronic properties, synthetic analogues that mimic these systems draw interest with potential applications in solar cells, sensors, or photodynamic therapy. 2−5 One of the systems that inspired a myriad of work is the self-assembled supramolecular structures of bacteriochlorophylls (BChls), which form highly efficient antennae to capture light without a protein scaffold to hold the structures together. 6−8 Replacing the central magnesium atom in BChl c with zinc has been a rewarding synthetic strategy to increase the stability of the metal complexes. 3 Moreover, by chemical modifications of the structures, their self-assembly processes and morphologies can be controlled. 9 Like their BChl counterparts, zinc-chlorins have rich intermolecular interaction types: π−π stacking, H-bonding, and metal−ligand. 10 An accurate description of intermolecular interactions is crucial for understanding the forces under which the supramolecular aggregates form. 11 Such descriptions also provide the necessary benchmark for improving the available force fields. 12 However, the studies on the intermolecular interactions governing the dimer formation of these zincporphyrins are scarce. Previously, Kuznetsov 13 analyzed groundstate geometries of Zn-porphyrin dimers and found good agreement with the experimental crystal structure. To the best of our knowledge, the relative importance of different types of interactions and the effect of the substitutions on the interaction energy (IE) strength are not known.
Here, for the first time, we present a systematic study of the IE in zinc-porphyrin complexes with increasing structural complexity so that the effect of substitutions on the dimer IEs is quantified. We used semiempirical, density functional, and symmetry-adapted perturbation methods to find a good balance of cost and accuracy. Moreover, the types of intermolecular interactions are evaluated using energy decomposition analysis based on the symmetry-adapted perturbation theory (SAPT).
A set of compounds starting from the simplest free-base porphyrin to more complex substituted Zn-chlorins are studied ( Figure 1). Although the focus is on the metal complexed structures, the free-base porphyrin (0) is included as a reference. Molecule 1 is a highly symmetrical Zn(II)-porphyrin for which an experimental crystal structure is available. 14 In molecule 2 (Zn-chlorin), the 17−18 bond is saturated, and hence, there are 20 π electrons. Molecule 3, similar to BChl c, has a hydroxy group attached to the 3 1 position, and the 13th, 14th, and 15th carbons are a part of a cyclopentenone. Molecule 4 has all of the methyl and ethyl group substitutions seen in BChl. This setup makes it possible to distinguish various groups' chemical and sterical contributions to the dimer IEs.

■ METHODS
For van der Waals complexes, the geometry relaxation upon complexation is small, changes in the interatomic distances are usually less than 10 −2 Å. Hence, the dimer IEs can be calculated with frozen monomer geometries. However, the amount of geometrical relaxation changes from one method to the other among the methods used here. For example, for some of the semiempirical methods, monomer interatomic distances change as much as 10 −1 Å upon dimer formation. In addition, unlike van der Waals complexes, the dimers studied have metal−ligand interactions with distances <3 Å. Therefore, both the binding energy (BE) and IE are calculated.
The BE, which is the total dimer energy with respect to the optimized monomer asymptote is The IEs for the dimers are calculated as follows Here, the monomer energies for their geometry in the dimer are subtracted from the total dimer energy. The basis set superposition error (BSSE) correction is obtained by the counterpoise method of Boys and Bernardi. 15 The performance of the density functional theory (DFT) methods is investigated for a number of functionals. The longrange corrected functional ωB97M-V 16 performed better than many others for a large number of data points. 16 Moreover, together with the def2-TZVP basis set, 17 it described the chemistry of transition metal complexes successfully. 18 For molecule 1, the results from the ωB97M-V functional is compared with the results from other density functionals from different rungs. These are B97-D3 19 from rung 1, TPSS-D3(BJ) 20,21 and B97M-rV 22 from rung 2, M06-2X 23 with and without -D3 dispersion, MN15, 24 ωB97X-D, 25 ωB97X-D3, 26 and ωB97X-V 27 from rung 3.
In addition, as low-cost alternatives to the methods mentioned above, the three corrections composite methods B97-3c 28 and PBEh-3c 29 and semiempirical methods DFTB3-D4 30,31 and GFN2-xTB 32 are investigated. B97-3c is a low-cost variant of B97-D that provides dispersion and avoids basis set errors with an optimized triple zeta basis (mTZVP). PBEh-3c yields geometries that are close to MP2/TZ accuracies with a fraction of the computational expense. DFTB3-D4 is a tightbinding approach that is also suitable for solid-state structure calculations. GFN2-xTB is also a semiempirical tight-binding method with multipole electrostatics and dispersion, which shows promise in the optimization of large aggregates and complexes of transition metals.
An essential technique for analyzing the noncovalent interactions is the SAPT. 33 The SAPT IE is a collection of the electrostatic, exchange, induction, and dispersion components. Within the SAPT formulation, under the assumption that intermolecular forces are much weaker than the intramolecular covalent interactions, the distortions in the geometry of a molecule in the presence of another one could be assumed negligible. Then, the IE is calculated as a perturbation in the monomer basis. If all of the higher-order terms in the perturbation expansion are included, it is possible to achieve accuracy close to that of CCSD(T). However, the standard level expansion (SAPT2) scales similarly to CCSD(T), ≈N 7 . So instead, here, the SAPT(DFT) 34 method together with the density-fitting approximation with frozen core orbitals 35 is used (the grac shift: 0.2498). Unfortunately, the def2-TZVP basis could be used only for molecules 0−2 due to the large basis set size and memory limitations (5478 basis functions for molecule 3).
The DFT calculations are performed as implemented in Q-Chem 5. To benchmark the accuracy of the computational methods, we first focus on the available Zn(II) porphyrin crystal structures. 14 Two crystal structures with P1̅ and P2 1 /c symmetries are available, both of which feature parallel displaced molecular packing in the form of dimer and trimers due to strong π−π stacking interactions. The monomers are relatively close-packed, with the distances between the planes being about 3.2 Å, which are shorter than what was reported for dimers of free-base porphyrin (3.4 Å), 45 and this perhaps is a result of the interaction of the metal atom with the nitrogen of the pyrrole ring in the neighboring monomer. This coordinative interaction pulls the metal off-center and closer to the other monomer, doming the N−Zn−N angles to 172−174°. The highly symmetric gas-phase structure of Zn-porphyrin (D 4h ) is no longer in the crystal due to intermolecular interactions. For example, there is an asymmetry in the intramolecular Zn−N bond distances.
First, we optimized the dimer structure extracted from the P2 1 /c crystal 14 by the ωB97M-V 16 method using the def2 basis family. 17 To make sure that the BSSE does not lead to an artificially high IE, we performed the dimer geometry optimization with def2-SVP, -TZVP, and -QZVP bases. 17 The difference in the IEs for triple and quadrupole zeta bases is small (Table 1). This finding confirms the results of earlier studies, which showed that the def2-TZVP basis is reliable without the counterpoise correction. 18,46 The BSSE in the IE for the dimer optimized at the ωB97M-V/ def2-TZVP level of theory is also calculated ( Table 2). As expected, the larger the basis set, the smaller the error. The BSSE with def2-TZVP is approximately 2 kcal/mol. Should the BSSE be assumed to be similar for both the BE and IE calculations, the difference between the IE and BE is the geometry relaxation energy due to the geometrical deformation of the monomers in the dimer. Hence, the relaxation energy is approximately 3 kcal/ mol compared to 2 kcal/mol BSSE. Therefore, ignoring the monomer geometry relaxation upon complex formation may not be justifiable if high accuracy is necessary.
Overall, the computed geometries show a good agreement with the crystal structures ( Table 1). Note that only the shortest Zn−N′ distances and the smallest N−Zn−N angles are reported in the tables. Similar to the crystal structure, the monomer is no longer symmetrical; there is a range of bond lengths and angles for the molecules in the dimer. For example, the intramolecular Zn−N distances vary in the range of 2.03−2.05 Å. In addition, there is a sufficient overlap of the optimized geometry for the dimer with the initial geometry extracted from the crystal. The high computational cost of using a quadruple-zeta basis is not justified since it does not significantly improve the geometry or the IE.
The geometry and BE of the dimer (molecule 1) obtained using ωB97M-V/def2-TZVP is compared with those obtained using other methods in Figure 2 and Table 3. All of the DFT results are approximately within 3 kcal/mol of each other. Previously, Kuznetsov studied the Zn-porphyrin dimer geometry using ωB97XD (Grimme's D2 dispersion model) with the def2-TZVP basis set and reported a BE in the same range (−27.3 kcal/mol) as that for molecule 1. 13 The approximate methods provide good accuracy while reducing the computational cost significantly ( Table 3). The "-3c" methods provide a good geometry, whereas the GFN2-xTB optimization brings the monomers too close together, reflected as the over-doming in the N−Zn−N angle, resulting in a larger IE. One notable feature in DFTB-D4 is the intramolecular Zn−N distances being too far. Overall, the parameterized methods provide good accuracy for geometry optimization. Most of the IEs lie between −25 and −30 kcal/ mol, the GFN2-xTB being the most significant deviation with a value of −35 kcal/mol. As a result of this comparison, to study the rest of the molecules, the following methods are chosen: B97M-rV, B97-D3, B97-3c, GFN2-xTB, and DFTB-D4.
For molecule 2, three parallel-stacked initial structures, which were defined according to the alignment of the saturated bond (i.e., θ = 0°when the saturated bond coincides in the two dimers) were used. The optimized geometries are close to these initial structures, but with parallel displacement similar to that in molecule 1 (Figure 3).
The optimized geometries at various levels of theory have similar BEs (Table 4). DFTB-D4 and GFN2-xTB predictions for the intermolecular Zn−Zn′ and N−Zn′ distances are slightly shorter, and the N−Zn−N angles have more doming. The GFN2-xTB BE is about 5 kcal/mol larger than the B97M-rV value, similar to molecule 1. Overall, the parameterized methods compare remarkably well with the B97M-rV/def2-TZVP level of calculations with a much lower cost. The similarity of the BE for rotated configurations of 2 indicates that a "sombrero hat" type of a potential energy surface is at play, as suggested previously by Muck-Lichtenfeld and Grimme for free-base porphyrin. 46 Molecule 3 includes both a hydroxy group and a ketone group; hence, in addition to π−π stacking, H-bonding and Zn− O interactions play a role in dimer formation. A total of six dimer geometries were used as the initial structures ( Figure 4). These were the four parallel-stacked dimers, rotated θ = 0, 90, 180, and 270°degrees around the Zn−Zn′ axis, and two paralleldisplaced geometries (initially displaced by about 6 Å) so that either hydroxy or keto oxygen on one monomer can interact with the metal atom on the other monomer. The 180°dimer with the double hydrogen bond has the strongest binding across the methods.
Again, optimized geometries of these dimers formed by rotation are all slightly sideways displaced, similar to molecules 1 The Journal of Physical Chemistry A pubs.acs.org/JPCA Article and 2. All of the final geometries are closer to the initial structure, indicating that the basin for the minima on the potential energy surface is wide. Consistently across the methods, the higher BE is for the dimer in which the θ = 180°rotation allows for two     Table 5. The structure of the θ = 180°H-bonded dimer with two hydrogen bonds is similar to the structure crystalized by Barkigia et al. 47 They observed that the dimers of Zn pheoporphyrin d orient to form nonlinear hydrogen bonds, similar to what we observed for molecule 3. (155°for the dimer optimized at the B97M-rV/TZVP level of theory). In addition, Barkigia et al. found that Zn−O distances in the complexes are larger than 2.1 Å, which is also the case for the optimized dimer geometries in Table 5.
The displaced dimers have weaker π−π interaction and no hydrogen bonding, but they are stabilized by Zn−O metal− ligand interaction. As a result, most of the optimized geometries are close to the starting geometries except for the B97-3c optimization for the displaced geometry with θ = 180°rotation, where the monomers slip back into the close pack form where the Zn−Zn′ distance is 3.42 Å. Surprisingly, although the geometries for displaced dimers optimized with the DFTB-D4 and GFN2-xTB methods are similar to those obtained using the B97M-rV method, the BEs are very large, almost the same as the BE for the face-to-face dimers. This could be attributed to an overestimate of the metal−ligand interactions; in the geometries optimized with tight binding methods, the Zn−O and Zn−N′ bonds are shorter and N−Zn−N angles are smaller (doming).
Similar to the dimers of molecules 1 and 2, the strongest interaction holding the molecule 3 dimers together is the dispersion, and consequently, the potential energy surfaces have wide minima. Therefore, any additional interaction such as metal−ligand or hydrogen bond and steric effects can modify the packing patterns, confirming the earlier findings of Muck-Lichtenfeld and Grimme for free-base porphyrin. 46 These authors noted that the IE differences between different motifs stem from the electrostatics. The same is true for the dimers of molecule 3. A previous work on the dimers of Chlorophyll a also had a similar conclusion. 10 IEs Depend on Substitutions. To quantify the effect of substitutions on the IE, the dimers corresponding to the lowest BE conformers of molecules 1−3 were analyzed using the SAPT(DFT) method. In particular, we used the 90 and 180°r otated dimers for molecules 2 and 3, respectively. In addition, we included two more dimers: (1) the free-base porphyrin dimer  Table 6. The decomposition of the IE into the components reveals that metal complexation and polar and R-group substitutions increase the IE ( Figure 5).
The BSSE-corrected IEs obtained using the B97M-rV/def2-TZVP level of theory are remarkably close to the SAPT(DFT)/ def2-TZVP values for molecules 1 and 2. On the other hand, although the def2-SVP basis underestimates the IE, the overall trends are more or less captured ( Figure 6). A comparison of the SAPT results with these two bases reveals that the major error in def2-SVP stems from the inadequate description of the dispersion interaction (Table 7).
For free-base porphyrin, Muck-Lichtenfeld and Grimme reported an IE value of −25 ± 3 kcal/mol. The values we obtained with def2-TZVP are similar, namely, −22 and −24 kcal/mol using DFT and SAPT methods, respectively.
Metal binding increases the dimer IEs by about −6 to −8 kcal/mol. We further analyzed the effect of metal binding in terms of the SAPT(DFT)/def2-TZVP components for molecules 1 and 2 (Table 7). Zn binding to the free-base porphyrin increases the electrostatics, dispersion, induction, and exchange interactions by −10.7, −8.5, −2.1, and 15.6 kcal/mol, respectively, resulting in an overall stabilization of dimers by approximately −6 kcal/mol. Furthermore, introducing a keto group and a hydroxy group further stabilizes the dimers by about −20 kcal/mol.

■ CONCLUSIONS
Zinc binding in free-base porphyrin stabilizes the dimer significantly. Introducing the hydroxy and keto groups to the Zn-chlorin ring almost doubles the IE compared to that of the free-base porphyrin. Interaction surfaces for all parallel stacked dimers have mostly flat minima, where the rotation does not significantly change the energies. These flat interaction surfaces are perhaps the reason behind rich packing patterns observed in experimentally produced aggregates. All complexes are dis-    The Journal of Physical Chemistry A pubs.acs.org/JPCA Article persion-bound, without which the BE would be positive. Although the dispersion is the strongest force, due to the flat PES, mostly electrostatics determines the geometries. Most methods we investigated provide reasonably good geometries and capture the IE for the face-to-face dimers, consistently predicting the doubly hydrogen-bonded dimer as the lowest energy configuration. However, GFN2-xTB and DFTB-D4 predict the IE for displaced dimers to be as high as that for the face-to-face dimers, whereas in the DFT methods, the IE is smaller for the displaced dimers, as expected from the diminishing dispersion interactions. On the other hand, B97M-rV performs comparably with SAPT(DFT) at the same basis set level. In addition, the analysis showed that ignoring the geometry relaxation causes a more significant error than the BSSE.
Remarkable charge and energy transport properties emerge only for supramolecular aggregates. Therefore, understanding the molecular interactions at the aggregate scale for various packing patterns is underway in our laboratory.